Accompanying Textbook:

The simplest data set encountered in functional data analysis is a sample of the form: \[ \begin{align} x_i(t_j),\quad t_j\in[a,b],\quad i=1,\dots,n \quad j=1\dots,J. \end{align} \] The theoretical objects we wish to study are smooth curves: \[ \begin{align} X_i(t),\quad t\in[a,b],\quad i=1,\dots,n \quad j=1\dots,J. \end{align} \]

Basis expansions:

Typically, the first step in working with functional data is to express them by means of a basis expansion \[ X_i(t)\approx\sum_{m=1}^M c_{im} B_m(t),\quad t\in[a,b], \] where \(B_1(t),\dots,B_M(t)\) are some standard collection of basis functions like:

• splines
• wavelets
• sine and cosine functions
• etc.

B-spline basis functions \(B_1(t),\dots,B_M(t)\):

bspl_bf <- create.bspline.basis(rangeval=c(0,1), norder=3, 
                                     breaks=seq(0,1,len=7))
plot(bspl_bf, main="B-spline Basis Functions", xlab="[a,b]=[0,1]")

Fourier basis functions \(B_1(t),\dots,B_M(t)\):

fourier_bf <- create.fourier.basis(rangeval=c(0,1), nbasis=5)
plot(fourier_bf, main="Fourier Basis Functions", xlab="[a,b]=[0,1]")

Example: Berkeley growth study raw data

matplot(x=growth$age, y=growth$hgtf[,1:5], type="p", lty=1, pch=21, 
        xlab="age", ylab="height (cm)", cex.lab=1.2,col=1:5,bg=1:5,
        main="5 Girls in Berkeley Growth Study")

Example: Berkeley growth study pre-processed data

SmObj <- smooth.basisPar(growth$age,y=growth$hgtf[,1:5],lam=0.1)
plot(SmObj$fd, xlab="age", ylab="height (cm)", cex.lab=1.2,
     main="5 Girls in Berkeley Growth Study", lty=1)

## [1] "done"

Example: Berkeley growth study 1st derivative

plot(deriv(SmObj$fd, 1), xlab="age", 
     ylab="growth rate (cm / year)", cex.lab=1.2,
     main="5 Girls in Berkeley Growth Study (1st derivative)", lty=1)

## [1] "done"

Example: Berkeley growth study 2nd derivative

plot(deriv(SmObj$fd, 2), xlab="age", cex.lab=1.2,
     ylab="growth acceleration (cm / year^2)",
     main="5 Girls in Berkeley Growth Study (2nd derivative)", lty=1)

## [1] "done"

Pointwise mean: \[ \bar{X}_n(t)=\frac{1}{n}\sum_{i=1}^n X_i(t) \] Pointwise standard deviation: \[ S_n(t)=\sqrt{\frac{1}{n-1}\sum_{i=1}^n\Big(X_i(t)-\bar{X}_n(t)\Big)^2} \]

Pointwise mean and standard deviation:

BM.Mean <- rowMeans(BM.mat)
BM.SD   <- apply(BM.mat, 1, sd)

Pointwise covariance function: \[ \hat{c}_n(t,s)=\frac{1}{n-1}\sum_{i=1}^n\Big(X_i(t)-\bar{X}(t)\Big)\Big(X_i(s)-\bar{X}(s)\Big) \]

• The sample covariance function is extensively used in Functional Data Analysis (FDA).
• The interpretation of the values of \(\hat{c}_n(t,s)\) is the same as for the usual variance-covariance matrix.

Pointwise covariance function:

BM.Cov <- var(BM.mat)

• Idea: Use the Estimated Functional Principal Components (EFPC's) \(\hat{v}_j\) as basis functions for \(X_i\): \[ X_i(t)\approx\bar{X}_n(t) + \sum_{j=1}^p\hat{\xi}_{ij}\hat{v}_j(t) \]

• Estimated scores: \(\hat{\xi}_{ij}=\int_a^b (X_i(t)-\bar{X}_n(t))\hat{v}_j(t)dt\)
• EFPC's \(\hat{v}_j\) are orthonormal, i.e. \[\textstyle\int_a^b\hat{v}_j(t)\hat{v}_k(t)dt=\begin{cases}1, &j=k\\0,&j\neq k.\end{cases}\]

• Best basis property

Eigendecomposition of \(\hat{c}_n(t,s)\): \[ \begin{align} \hat{c}_n(t,s)&\approx\sum_{j=1}^p\hat{\lambda}_j\hat{v}_j(t)\hat{v}_j(s), \end{align} \] where \(\hat{\lambda}_1>\dots >\hat{\lambda}_p\) denote the estimated eigenvalues.

\[ \begin{align} X_i(t)&\approx\bar{X}_n(t) + \sum_{j=1}^p\hat{\xi}_{ij}\hat{v}_j(t) \end{align} \]

• Bank of America (BOA) stock values \(P_i(t_j)\)
  • \(i=1,\dots,n=2511\) denotes the trading days
    (from April 9th, 1997 to April 2nd, 2007).
  • \(t_j\in[0,6.5]\) denotes intra-day trading time
    (6.5 hours / per trading day).
  • Minute-wise measurements: \(j=1,\dots,J=390\).

We study the daily cumulative log-return functions: \[R_i(t):=\log(P_i(t))-\log(P_i(0))\approx\frac{P_i(t)-P_i(0)}{P_i(0)}\]

# Plot functional data
plot(log_BOA_fd, xlab="Trading Hours",ylab="",lwd=1, col=gray(.5), 
     main="Cumulative Log-Return Functions")
lines(log_BOA_fd[1:10],lwd=1.5, lty=1)

• A function \(f\) is said to be square integrable, i.e., \(f\in L^2([a,b])\) if \[ \textstyle\int_a^b \big(f(t)\big)^2 dt<\infty. \]
• Square integrable functions form a vector space, i.e., for \(c,d\in\mathbb{R}\): \[g,f\in L^2([a,b])\quad\Rightarrow\quad cf+dg\in L^2([a,b]),\] where addition is pointwise, i.e.
  \((cf+dg)(t)=cf(t)+dg(t)\) for almost all \(t\in[0,1]\).
  
  (We will ignore measure-theoretic considerations.)

What makes the space \(L^2([a,b])\) so convenient is its structure:

• Inner product of two functions is given by \[ \textstyle\langle f,g\rangle=\int_a^bf(t)g(t)dt \]

• Two functions \(f\) and \(g\) are orthogonal if \[\langle f,g\rangle=0\]

• The norm \(||f||=\sqrt{\langle f,f\rangle}\), gives us a notion for the distance between two functions: \[d(f,g)=||f-g||\]

As we have already seen, basis expansions play an important role in methodology for functional data.

• We say that a set of functions \(\{e_1 ,e_2,e_3,\dots\}\) is a basis in \(L^2([a,b])\) if every function \(f\in L^2([a,b])\) admits a unique expansion \[ \textstyle f(t)=\sum_{j=1}^\infty a_j e_j(t) \]
• We say that \(\{e_1 ,e_2,e_3,\dots\}\) is an orthonormal basis if, in addition \[\langle e_j , e_k\rangle=\begin{cases}1, &j=k\\0&j\neq k.\end{cases} \]

Let \(X\) denote a random function defined on a probability space, say \(\Omega\).

• We assume that all realizations \(X(\omega)\), \(\omega\in\Omega\), are elements of the space \(L^2([a,b])\) of square integrable functions, i.e
  
  \(||X(\omega)||=\sqrt{\int_a^b \big(X(\omega)(t)\big)^2dt} < \infty\) for all \(\omega\in\Omega\).

• \(||X||\in\mathbb{R}\) is thus a real random variable.

• If \(E(||X||^2)<\infty\), we say that the random function \(X\) is square integrable. Caution: Integration is here with respect to \(\omega\), not \(t\). It might be more pedagogical to say that the random function \(X\) has a finite second moment.

• Mean function: \[ \mu(t)=E(X(t)) \]
• Covariance function: \[ c(t,s)=E\Big((X(t)-\mu(t))(X(s)-\mu(s))\Big) \]
• The sample mean function \(\bar{X}_n(t)\) and the sample covariance function \(\hat{c}_n(t,s)\) are viewed as estimators of the population parameters \(\mu(t)\) and \(c(t,s)\).

The population covariance function leads to Functional Principal Component Analysis (FPCA), which allows us to represent a square integrable random function \(X\) as: \[ \textstyle X(t)=\mu(t)+\sum_{j=1}^\infty\xi_j v_j(t) \]

• The eigenfunctions \(v_j\) are the solutions of the eigen-equation \(\int_a^bc(t,s)v_j(s)ds=\lambda_jv_j(t)\).

• \(\lambda_1\geq\lambda_2\geq\dots\) denote the eigenvalues.

• The random variables \(\xi_j\) are called the scores \[ \textstyle \xi_j=\langle X-\mu,v_j\rangle=\int_a^b(X(t)-\mu(t))v_j(t)dt \]

FPCA allows us to represent a square integrable random function \(X\) as ( Karhunen-Loéve expansion): \[ \textstyle X(t)=\mu(t)+\sum_{j=1}^\infty\xi_j v_j(t) \]

It can be shown that: \[ E(\xi_j)=0,\quad V(\xi_j)=\lambda_j,\quad \operatorname{Cov}(\xi_j,\xi_k)=0,\,j\neq k, \] and \(E\int_a^b(X(t)-\mu(t))^2dt=E||X-\mu||^2=\sum_{j=1}^\infty\lambda_j\).

That is, \(\lambda_j\) is the variance of the random function \(X\) in the principal direction \(v_j\) and the sum of these variances is the total variance of \(X\).

Scalar–on–function regression: \[Y_i=\int_a^b\beta(s)X_i(s)ds+\varepsilon_i\] Function–on-scalar regression: \[Y_i(t)=\sum_{k=1}^p\beta_k(t)x_{ik}+\varepsilon_i(t)\] Function–on-function regression: \[Y_i(t)=\int_a^b\beta(t,s)X_{i}(s)ds+\varepsilon_i(t)\]

Nonparametric scalar–on–function regression: \[Y_i=m(X_i)+\varepsilon_i,\] where \(Y_i\in\mathbb{R}\) and \(X_i\in L^2([a,b])\).

Nonparametric function–on–function regression: \[Y_i=m(X_i)+\varepsilon_i,\] where \(Y_i,\,X_i\in L^2([a,b])\).

The above models are just prototypes illustrating the general idea. The main point is that the functional regression models are infinite dimensional objects which must be estimated from a finite sample.

Remember: The usual regression model \[ \textstyle Y_i=\sum_{r=1}^p\beta_rX_r+\varepsilon \] Normal equations and solutions (population & sample): \[ \textstyle \mathbf{C}_{XY}=\mathbf{C}_{XX}\boldsymbol{\beta}\quad\Rightarrow\quad \boldsymbol{\beta}=\mathbf{C}_{XX}^{-1}\mathbf{C}_{XY} \] \[ \textstyle \widehat{\mathbf{C}}_{XY}=\widehat{\mathbf{C}}_{XX}\widehat{\boldsymbol{\beta}}_n\quad\Rightarrow\quad \widehat{\boldsymbol{\beta}}_n=\widehat{\mathbf{C}}_{XX}^{-1}\widehat{\mathbf{C}}_{XY} \] where
\([\mathbf{C}_{XX}]_{r,\ell}=E(X_rX_\ell)\),
\([\mathbf{C}_{XY}]_r=E(X_rY)\),
\(\big[\widehat{\mathbf{C}}_{XX}\big]_{r,\ell}=\frac{1}{n}\sum_{i=1}^n X_{ir}X_{i\ell}\), and
\(\big[\widehat{\mathbf{C}}_{XY}\big]_r=\frac{1}{n}\sum_{i=1}^n X_{ir}Y_{i}\).

Normal equations and solutions (population & sample): \[ \begin{align} \textstyle \mathbf{C}_{XY}=\mathbf{C}_{XX}\boldsymbol{\beta}\quad&\Rightarrow\quad \boldsymbol{\beta}=\mathbf{C}_{XX}^{-1}\mathbf{C}_{XY}\\ \textstyle \widehat{\mathbf{C}}_{XY}=\widehat{\mathbf{C}}_{XX}\widehat{\boldsymbol{\beta}}_n\quad&\Rightarrow\quad \widehat{\boldsymbol{\beta}}_n=\widehat{\mathbf{C}}_{XX}^{-1}\widehat{\mathbf{C}}_{XY} \end{align} \]

To obtain \(\boldsymbol{\beta}\) and \(\widehat{\boldsymbol{\beta}}_n\), we need invertibility of the \(p\times p\) matrices \(\mathbf{C}_{XX}\) and \(\widehat{\mathbf{C}}_{XX}\).

Functional regression \[ \textstyle Y_i=\int_a^b\beta(t)X_i(t)+\varepsilon_i \] Normal equations (population & sample): \[ \textstyle c_{XY}(t)=\int_a^b c_{XX}(t,s)\beta(s)ds \] \[ \textstyle \hat{c}_{XY}(t)=\int_a^b \hat{c}_{XX}(t,s)\beta(s)ds \] where
\(c_{XX}(t,s)=E(X(t)X(s))\),
\(c_{XY}(t)=E(X(t)Y)\), \(\hat{c}_{XX}(t,s)=\frac{1}{n}\sum_{i=1}^n(X_i(t)X_i(s))\), and
\(\hat{c}_{XY}(t)=\frac{1}{n}\sum_{i=1}^n(X_i(t)Y_i)\).

Functional regression \[ \textstyle Y_i=\int_a^b\beta(t)X_i(t)+\varepsilon_i \] Normal equations (sample): \[ \textstyle \hat{c}_{XY}(t)=\int_a^b \hat{c}_{XX}(t,s)\beta(s)ds \]

Problem: - Ill-posed inversion problem - Eigenvalues \(\approx 0\)

Classical solutions: - Ridge regression - Principal component regression

Let us consider the scalar-on-function model \[ \textstyle Y_i=\alpha + \int_a^b \beta(s)X_i(s)ds+\varepsilon_i,\quad i=1,\dots,n \]

Simplest approach: expand \(\beta\) using deterministic basis functions (splines, fourier, etc) \[ \textstyle \beta(t)=\sum_{k=1}^Kc_kB_k(t) \]

Rewrite the functional regression model into an approximative multivariate reggression problem \[ \begin{align} \int_a^b\beta(s)X_i(s)ds &\approx\sum_{k=1}^Kc_k\int_a^bB_k(s)X_i(s)ds\\ &=\sum_{k=1}^Kc_kx_{ik} \end{align} \]

Let us consider the scalar-on-function model